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Constraints on a planetary origin for the gap in GM 
Aurigae's protoplanetary disc 
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ABSTRACT 

The unusual spectral energy distribution (SED) of the classical T Tauri star GM 
Aurigae provides evidence for the presence of an inner disc hole extending to several 
au. Using a combination of hydrodynamical simulations and Monte Carlo radiative 
transport, we investigate whether the observed SED is consistent with the inner hole 
being created and maintained by an orbiting planet. We show that an ~ 2Afj up i te r 
planet, orbiting at 2.5 au in a disc with mass O.O47M and radius 300 au, provides a 
good match both to the SED and to CO observations which constrain the velocity field 
in the disc. A range of planet masses is allowed by current data, but could in principle 
be distinguished with further observations between 3 and ~ 20 microns. Future high 
precision astrometric instruments should also be able to detect the motion of the 
central star due to an orbiting Jupiter mass planet. We argue that the small number 
of T Tauri stars with SEDs resembling that of GM Aur is broadly consistent with the 
expected statistics of embedded migrating planets. 
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1 INTRODUCTION 

The detection of extrasolar planets and planetary systems 
(e.g., Mayor & Queloz 1995) has invigorated studies of plan- 
etary formation (Lissauer 1993; Boss 1998, 2000), migration 
(e.g., Lin & Papaloizou 1986; Lin et al. 1996; Trilling et al. 
1998; R. Nelson et al. 2000), and interaction with a pro- 
toplanetary disc (Lin & Papalizou 1979a, b; Artymowicz 
& Lubow 1994). Recent theoretical work has investigated 
mechanisms for halting planetary migration and explaining 
the statistical distribution of the orbital radii of extraso- 
lar planets (Armitage et al. 2002; Kuchner & Lecar 2002; 
Trilling, Lunine & Benz 2002) . Theoretical models for the in- 
teractions of (proto)planets with a protoplanetary disc show 
that large low density regions (gaps) can be created within 
the disc (Lin & Papaloizou 1979a, b; Artymowicz & Lubow 
1994). The detection of these gaps would provide strong sup- 
port for planetary formation theories, but current observa- 
tions do not have the spatial resolution to directly image 
the au sized gaps that are predicted. Future instrumenta- 
tion, such as ALMA, may allow direct detection of gaps 
(Wolf, Henning & Kley 2002), but for now we are restricted 
to spectral techniques for inferring the presence of disc gaps. 



Spectral energy distributions (SEDs) from protoplan- 
etary discs arise from the thermal reprocessing of starlight 
and dissipation of viscous accretion luminosity (Lynden-Bell 
& Pringle 1974; Adams, Lada & Shu 1987; Kenyon & Hart- 
mann 1987). The resulting SEDs display characteristic in- 
frared excesses that are sensitive to the size, shape, and mass 
of the disc. Disc gaps may manifest themselves in SEDs by 
removing material that would normally have contributed to 
the SED at a given wavelength: gaps close to the star re- 
sult in smaller than usual near-IR excesses, while gaps at 
large disc radii will cause distortions to the SED at longer 
wavelengths (Beckwith 1999). The classical T Tauri stars 
GM Aur and TW Hya have SEDs that indicate the inner 
regions of their discs have been cleared of material out to 
around 4 au from the central stars (Koerner, Sargent, & 
Beckwith 1993; Chiang & Goldreich 1999; Schneider et al. 
2002; Calvet et al. 2002). 

There are a number of possible mechanisms for clear- 
ling such gaps. Photoevaporation by ultraviolet radiation, 
either from the central star or from surrounding radiation, 
can produce a disc wind at radii where the sound speed 
exceeds the escape velocity (Clarke et al. 2001). The inner 
regions of the disc will then be cleared on a viscous timescale 
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once the wind mass loss rate is comparable to the accretion 
rate. This, however, generally requires an accretion rate two 
orders of magnitude lower than the lO~ 8 MQyr~ 1 accretion 
rate expected for GM Aur (Gullbring et al. 1998). Clarke et 
al. (2001) also assumed a constant ionizing flux. Matsuyama 
et al. (2002) show that the clearing of the inner disc is even 
less likely if the decrease in ionizing flux with decreasing ac- 
cretion rate is included in the calculation. A. Nelson et al. 
(2000) also suggest that viscous dissipation and the dissipa- 
tion of gravitational instabilities may increase the midplane 
temperature of the inner regions of the disc to above the 
grain destruction temperature. This would then manifest it- 
self as a reduction in the flux of radiation with wavelengths 
between a few and ten microns. Their simulation, however, 
assumed a disc luminosity (0.5Lq) an order of magnitude 
greater than that calculated by Gullbring et al. (1998) for 
GM Aur (0.071L Q ). 

It has been suggested (Calvet et al. 2002) that the inner 
regions of TW Hydra's disc could be cleared by tidal interac- 
tions with an embedded planet (Lin & Papalizou 1979a, b). 
In this paper we consider if the gap in GM Aur's disc could 
be formed in the same way. We use a dynamical model of an 
approximately Jupiter mass planet in a protoplanetary ac- 
cretion disc, with parameters appropriate for GM Aur (e.g., 
Koerner, Sargent & Beckwith 1993; Simon, Dutrey & Guil- 
loteau 2000; Schneider et al. 2002), to investigate if the re- 
sulting disc structure can reproduce the observed SED. We 
also consider other planet masses to test if the SED can be 
used to constrain the planet mass. Estimates for GM Aur's 
age vary from ~ 1.5 Myr (Beckwith et al. 1990; Siess, Fores- 
tini & Bertout 1999) to ~ 10 Myr (Simon & Prato 1995; 
Hartmann 2002). There has therefore been sufficient time for 
an embedded planet to have formed via core accretion (Lis- 
sauer 1993) or via gravitational collapse (Boss 1998, 2000). 
The observation of a much younger system (e.g., < 1 Myr) 
with an SED signature indicating the presence of a disc gap 
would, if the gap could be shown to be due to the presence 
of an embedded planet, suggest that either core accretion 
occurs on a much shorter timescale than expected or that 
planet formation can occur via an alternative mechanism 
such as gravitational collapse. In section 2 we describe and 
present the results of our dynamical simulations, in § 3 we 
compare model SEDs with the observed SED of GM Aur, 
and we summarize our results in § 4. 



2 DYNAMICAL SIMULATIONS 

2.1 Dynamical Model and Initial conditions 

The simulations presented here were performed using 
smoothed particle hydrodynamics (SPH), a Lagrangian hy- 
drodynamics code (e.g., Benz 1990; Monaghan 1992). We 
initially assume a planet with mass M p = 1.7Mj orbits a 
star, with mass M* = O.85M0, at a radius of r p — 2.5 
au. Both the central star and planet are modelled as point 
masses onto which gas particles can accrete if they approach 
to within the sink radius (Bate et al. 1995). The planet or- 
bits within a circumstellar disc of mass Mdisc = O.O47M0, 
modelled using 300000 SPH gas particles distributed such as 
to give a surface density profile of E oc r _1 (D'Alessio et al. 
1998). The disc extends from n n = 0.25 au to r out = 300 au. 



The disc temperature varies as T oc r -1 / 2 , as expected for 
discs heated by starlight (D'Alessio et al. 1998), and we as- 
sume hydrostatic equilibrium in the vertical direction, giving 
a central density profile of p oc r~ 2 ' 25 . For the SPH simu- 
lation, the disc is assumed to be locally isothermal and the 
above temperature profile is maintained throughout the sim- 
ulation. The parameters chosen are appropriate for GM Aur 
based on calculations of the stellar mass (Simon, Dutrey & 
Guilloteau 2000) and the disc mass and size (Schneider et 
al. 2002). 

The vertical scale height, H, of the disc is given by 
H 2 — c 2 r 3 /GM, where c s is the local sound speed, r is 
the radius, G is the gravitational constant, and M„ is the 
stellar mass. The temperature is normalised such as to give 
a scale height of 10.5 au at 100 au. The scale height at 2.5 
au is then 0.1 au, giving an H/r value of 0.04. This choice of 
scale height was based on an SED model using an analytic 
disc geometry and is similar to that obtained by Schneider 
et al.(2002). In using the above form for the scale height, we 
assume that the dust has the same scale height as the gas. 
Although this may not necessarily be the case, the need to 
use a dust scale height that has the same form as a standard 
equilibrium disc (e.g., Pringle 1981) would seem to make this 
a reasonable assumption. 

Using 300000 SPH particles the vertical profile of the 
disk can be modelled accurately within ~ 3 scale heights 
of the disc midplane. For higher altitudes, the density is 
calculated analytically using (e.g., Pringle 1981) 

p(*,r)=p c (r)exp(-z 2 /2H 2 ) (1) 

where z is the height above the disc midplane, and p c (r) 
is the central density a distance r from the central star, 
determined using the dynamical model. Very little mass is 
contained in these regions of the disc and is hence unimpor- 
tant for the dynamical calculation. It is, however, extremely 
important for calculating the radiative transfer through the 
disc (see D'Alessio et al. 1998). 

Self-gravity is included in the calculation and both the 
point masses and the gas use a tree to calculate gravitational 
forces (Benz et al. 1990). The inclusion of self-gravity means 
that the central star and planet are both free to move under 
the gravitational influence of the disc. A great saving in 
computational time is obtained by using individual, particle 
time-steps (Bate et al. 1995; Navarro & White 1993) which 
are limited by the Courant condition and a force condition 
(Benz et al. 1990). 

2.2 Planet-disc interactions 

A planet orbiting within a circumstellar disc will tidally in- 
teract with the disc. Disc material exterior to the planet 
will gain angular momentum from the planet and will move 
to larger radii, while material interior to the planet will 
lose angular momentum and move to smaller radii (Lin & 
Papaloizou 1979a, b; Artymowicz & Lubow 1994). Viscous 
processes within the disc will tend to act against the tidal 
forces and the formation of a gap will depend on the tidal 
forces overcoming these viscous processes (Lin & Papaloizou 
1979a). This depends on the relative magnitudes of the mass 
ratio of the planet and the star, q, and the Reynolds num- 
ber within the disc, TZ. If q > TiT 1 a gap will form, while if 
q 2 > VT 1 the disc will be truncated at the outer Lindblad 
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resonance. The Reynolds number is given by Q.r 2 /u where 
f2 is the angular frequency, r is the radius, and v is the tur- 
bulent viscosity in the disc generally described as an alpha 
viscosity (Shakura & Sunyaev 1973). 

For the initial simulation (M v = 1.7Mj, Mdi sc = 
O.O47M ) considered here q — 0.002, and the Reynolds num- 
ber, determined from the artificial viscosity included in the 
code (Murray 1996), is ~ 20000. Therefore q > Tl' 1 but 
q 2 < VT 1 . A gap should form with the gap edge at the 
outer Lindblad resonance located, for a planet orbiting at 
2.5 au, at a radius r — 4 au. However, for the mass ra- 
tio and Reynolds number in this initial simulation, the disc 
should not be strongly truncated at the outer Lindblad res- 
onance. The Lindblad resonance can be determined from 
f2(r) = wm/(l + m) where ui is the angluar frequency of 
the orbiting planet and, for the outer Lindblad resonance, 
m — 1 (Binney & Tremaine 1987; Lin & Paploizou 1979a). 
The angular frequency at the outer Lindblad resonance is 
therefore f2(r) = cu/2. The assumption of a planet orbiting 
at 2.5 au was based on SED calculations using an analytic 
disc geometry and the results of Schneider et al.(2002) which 
suggest that the best fit to the SED occurs for gap extend- 
ing to a radius of 4 au (i.e., the radius of the outer Lindblad 
resonance for a planet orbiting at 2.5 au). 

Figure 1 shows the azimuthally averaged midplane den- 
sity at the beginning (dashed line) and end (solid line) of the 
simulation. After 1820 years the planet orbiting at 2.5 au has 
produced a gap in the disc out to r ~ 4 au. The gap edge 
is relatively soft, consistent with the relative magnitudes of 
the mass ratio (q) and Reynolds number (TV). The softness 
of the gap edge is likely to have an impact on the radiation 
transfer through the disc. As will be considered in a later 
section, to produce a gap with a significantly harder edge 
would require a more massive companion or a significantly 
less viscous disc. The mass within 5 au has reduced from 
4 x 1O~ 4 M at t = to 4 x 10~ 5 M Q at t = 1820 years. The 
reduction in the number of SPH particles respresenting the 
gas within 5 au requires the use of equation 1 to determine 
the vertical density at distances greater than 1 scale height 
above the disc midplane. 

If the planet clears sufficient material from within its 
orbit, torques from material outside its orbit should cause 
it to eventually lose angular momentum and spiral inwards. 
However, the planet in this simulation is approximately four 
times as massive as the disc material contained within its 
orbit and hence we expect it to suffer minimal orbital mi- 
gration (Syer & Clarke 1995; R. Nelson et al. 2000). 

Accretion onto the central star and the orbiting planet 
continues throughout the simulation despite the formation 
of the gap, consistent with simulations performed by Arty- 
mowicz & Lubow (1996). The central star has increased 
from it's initial mass of 0.85M Q to a mass of 0.85019M Q 
after 1820 years, while the orbiting planet has increased 
from 1.7 x 10~ 3 M Q to 1.88 x 1O~ 3 M in the same time 
period. The accretion rate at the end of the simulation (av- 
eraged over the last 100 years) was ~5x lO~ 8 M0yr _1 and 
1 x lO _8 M0yr _1 for the central star and orbiting planet 
respectively. An accurate calculation of the accretion rate 
in our simulation is limited by the mass resolution, since 
each SPH particle has a mass of 1.6 x 1O~ 7 M , and by our 
simplistic accretion process. We do not include the rotating 
magnetosphere of the T Tauri star which is likely to trun- 
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Figure 1. Azimuthally averaged midplane disc density at the 
beginning (dashed line) of our simulations and after the 1.7Mj 
planet has completed 450 orbits at 2.5 au (solid line). 



cate the disc at the corotation radius and force the accretion 
flow to follow the stellar magnetic field lines (Ghosh & Lamb 
1978). The accretion rate from the model is however similar 
to that expected for GM Aur (Gullbring et al. 1998). That 
mass is flowing through the gap and accreting onto the cen- 
tral star and planet is also consistent with other simulations 
(Artymowicz & Lubow 1996; Lubow et al. 1999). 

The disc velocities also remain very close to Keplerian 
despite the presence of an orbiting planet that perturbs the 
potential. This is consistent with CO observations (Koerner, 
Sargent & Beckwith 1993; Dutrey et al. 1998) which suggest 
that GM Aur is surrounded by a Keplerian circumstellar disc 
inclined at ~ 50°. 



3 RADIATION TRANSFER MODELS 

We calculate the SED for our simulations using the Monte 
Carlo radiative equilibrium technique of Bjorkman & Wood 
(2001). The output of our radiative transfer simulation is the 
disc temperature structure and the emergent SED. The disc 
heating is assumed to be from starlight and the Monte Carlo 
"photons" are tracked within a two dimensional grid in r and 
9. Our radiation transfer code calculates disc temperature 
and does not assume a vertically isothermal structure as in 
the SPH simulation. The input stellar spectrum is a 4000K, 
loggr = 3.5 Kurucz model atmosphere (Kurucz 1994) and 
the stellar radius (R*) is taken to be 1.75 Rq. The circum- 
stellar dust opacity is taken to be that described in Wood 
et al. (2002a) which fits the HH 30 IRS SED and was also 
used to fit the GM Aur SED (Schneider et al. 2002). The 
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dust model has a shallower wavelength dependent opacity 
than "standard" ISM dust (e.g., Mathis, Rumpl, & Nord- 
sieck 1977) and a larger grain size as inferred for many other 
classical T Tauri stars (e.g., Beckwith et al. 1990; Beckwith 
& Sargent 1991). 

We find that the dearth of near-IR excess in GM Aur's 
SED cannot be explained by a dust composition different to 
that used here. Using an ISM mixture and a disc without 
a substantial gap still results in a near-IR excess not ob- 
served in GM Aur. Our results are also in agreement with 
the Calvet et al. (2002) modelling of TW Hya. In their model 
they do use different dust types in the inner and outer disc, 
but still require a very low density inner disc to reproduce 
the near-IR emission observed in TW Hya. Further, ISM 
grains yield too steep a slope to the sub-mm SED (see also 
Beckwith & Sargent 1991). We therefore conclude that the 
dearth of near-IR excess emission from GM Aur can only be 
explained by an evacuated (low density) inner region. 

Circumstellar discs absorb and reprocess stellar pho- 
tons in a layer which, depending on the disc optical depth, 
can lie four to five scaleheights above the disc midplane 
(D'Alessio et al. 1999). Our SPH simulations cannot re- 
solve this low density material, so to provide a density grid 
for our radiation transfer simulations we have extrapolated 
the density, using Equation (1), to large heights above the 
disk midplane. The density calculated from the SPH simula- 
tions was azimuthally averaged and then interpolated onto 
a two dimensional grid in r and 8. The grid has 100 cells for 
0.25au < r < 300au spaced proportional to r 3 and 400 cells 
linearly spaced for < 9 < tv. This gridding allows us to ad- 
equately resolve the SPH density for our SED calculations. 

3.1 GM Aur SED 

Figure 2 shows the GM Aur SED data. The dotted line 
shows the input model atmosphere, and the open squares are 
data points obtained using various instruments (e.g., Cohen 
& Kuhi 1979; Kenyon & Hartmann 1995; Weaver & Jones 
1992; Beckwith & Sargent 1991). The solid line shows the 
SED calculated using the azimuthally averaged SPH density 
as input to the Monte Carlo radiation transfer calculation. 
The SPH simulation was run for 450 orbits of the 1.7 Jupiter 
mass planet at 2.5 au (1820 years), giving sufficient time for 
the density profile to reach a steady state. The disc struc- 
ture is also approximately azimuthally symmetric since the 
viscous timescale is significantly greater than the planet's or- 
bital period. The inclination of the disc was taken to be 50° , 
consistent with the value of 54° ± 5° determined by Dutrey 
et al. (1998), and with the recent scattered light models of 
Schneider et al. (2002). The model SED provides a good fit 
to the observations, particularly the lack of near-IR excess 
and the level and slope of the long wavelength continuum. 
The dashed line in Figure 2 is the SED for an equivalent 
disc in which a gap is not present and in which the inner 
edge of the disc is taken to be at 7 stellar radii. Comparison 
between the solid line and dashed line clearly shows that the 
effect of the disc clearing is to reduce the excess emission at 
short wavelengths (2 — 10 microns) and slightly increase the 
excess at longer wavelengths (10 — 30 microns). 

This is consequence of the stellar flux being reprocessed 
at the disc edge. In the simulation with i? m i n = 7R* (dashed 
line in Figure 2) stellar photons are reprocessed in high tem- 




A( ( un) 

Figure 2. GM Aur's spectral energy distribution (squares), our 
model viewed at i = 50° (solid line) and input stellar spectrum 
(dotted line). The dashed line shows the SED calculated in the 
absence of a gap with the disc extending in to 7 stellar radii. The 
presence of a gap (solid line) clearly reduces the near infra-red 
flux and enhances the flux at slightly longer wavelengths. 

perature (> 1000 K) regions of the disc giving a near-IR 
excess. When a gap is present the photons are reprocessed 
further out in the disc where the temperature is lower, de- 
creasing the near infrared flux and increasing the flux at 
longer wavelengths. A fit to the data is only possible if the 
gap extends all the way to the inner disc boundary. If the 
planet has only had sufficient time to clear a ring of material, 
leaving the inner disc almost unchanged, this inner material 
should remain hot, resulting in 2 and 3 micron fluxes similar 
to that for the dashed line in Figure 2, and higher than that 
observed for GM Aur (open squares). 

3.2 Companion (planet) mass and viscosity 

Although we are able to reproduce the observed SED with a 
1.7 Jupiter mass planet orbiting at 2.5 au, the exact form of 
the density profile depends on the relationship between the 
mass ratio, q, and the disc Reynolds number, 1Z. A weakness 
of SPH is the difficulty in determining exactly the viscosity 
and hence Reynolds number. In light of this, we have per- 
formed a number of additional simulations considering sub- 
stellar companions (planets) of various masses, each of which 
was run for at least 250 orbits of the companion around the 
central star (> 1000 years). In these simulations the disc 
was modelled between 0.25 and 100 au using 150000 SPH 
particles with a surface density profile such that the total 
disc mass would be O.O47M0 were the disc extended to 300 
au. Using 150000 SPH particles to simulate the inner 100 au 
of the disc should give approximately the same resolution as 
the 300000 particles used to simulate the disc out to 300 au. 
The Reynolds number in the simulations using 150000 SPH 
particles should therefore be very similar to the Reynolds 
number in a 300000 particle simulation (Murray 1996). The 
companion masses considered for these additional simula- 
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tions were 0.085, 21, and 43Mj. In simulating companions 
with relatively large masses (21 and 43Mj) we are not neces- 
sarily suggesting that one would expect to observe substellar 
companions with such masses, but simply illustrating how 
the SED can vary with mass. Strictly speaking, such mas- 
sive objects would generally be regarded as brown dwarfs, 
and we therefore refer to the 21 and 43Mj companions as 
substellar objects rather than planets. 

The process of star formation is known to preferen- 
tially produce binary stellar systems rather than single stars 
(Duquennoy & Mayor 1991). On the main sequence, how- 
ever, the mass ratio of close binaries is rarely so extreme 
as to pair a Solar mass star with a brown dwarf (Marcy & 
Butler 2000; Halbwachs et al. 2000). If this trend also exists 
during the pre-main-sequence phase, one would expect disc- 
embedded planetary mass objects to be more common than 
embedded close brown dwarfs. The formation mechanism 
of close binaries is, however, unclear with some simulations 
(Bate, Bonnell & Bromm 2002) suggesting substantial evo- 
lution of binary orbital parameters while a disc is present. 
The possibility of observing substellar mass companions at 
this early stage could therefore be an important test of dif- 
ferent binary formation models. 

The Reynolds number, 1Z, in the disc is calculated to be 
~ 20000. The 0.085 and 1.7 Jupiter mass planets therefore 
have q > TIT 1 but q 2 < TZ^ 1 while the 21 and 43 Jupiter 
mass substellar objects have q > VT 1 and q 2 > TIT 1 . The 
two planets (0.085 and 1.7Mj) should therefore produce 
gaps that are not strongly truncated at the outer Lindblad 
resonance while the two substellar objects (21 and 43Mj) 
should truncate the disc strongly at the outer Lindblad res- 
onance (Lin & Papaloizou 1979a; Syer & Clarke 1995). For 
a Reynolds number larger than used here, a planetary mass 
companion (< 10Mj) would satisfy easily the latter condi- 
tion (q > TIT 1 and q 2 > IZ^ 1 ). The results we obtain for the 
companions referred to as substellar objects would therefore 
also be applicable to planetary mass objects in a disc with 
a Reynolds number larger than in the simulation presented 
here. 

Figure 3 shows the azimuthally averaged midplane den- 
sities for the three additional masses together with the az- 
imuthally averaged midplane density profile for the 1.7Mj 
planet, between radii of 0.25 and 25 au. As expected the 
0.085 and 1.7 Jupiter mass planets produce gaps with soft 
edges while the 21 and 43 Jupiter mass objects truncate the 
disc strongly at the outer Lindblad resonance. The position 
of the gap edge varies slightly for the different masses, but 
in all cases is close to the outer Lindblad resonance at ~ 4 
au. 

Figure 4 shows the SEDs calculated using the az- 
imuthally averaged density profiles for the four masses as 
input to the radiation transfer calculation. Since the sim- 
ulations considering masses of 0.085, 21, and 43Mj were 
only performed out to 100 au, the density profiles were ex- 
trapolated to 300 au. The solid and dotted lines are for the 
1.7M.J and 0.085Mj planets, while the dashed and dash- 
dot lines are for the 43Mj and 21Mj substellar objects. 
The SEDs for the 1.7 and 0.085Mj planets are very simi- 
lar, as are the SEDs for the two substellar objects (21 and 
43Mj), which are, in fact, very difficult to separate. Figure 
5 shows a detailed view of the computed SEDs between 5 
and 200 microns with the line styles the same as in Figure 
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Figure 3. Azimuthally averaged midplane density profiles for 
substellar objects (planets) with masses of 0.085 (dotted line), 
1.7 (solid line), 21 (dash-dot line), and 50 Mj (dashed line). The 
two higher mass objects truncate the disc more strongly at the 
outer Lindblad resonance than the two lower mass planets. 



4. Figure 5 also shows ±25 % error bars on the data points 
(open squares). This choice of error was based on worst case 
1 a error estimates (Cohen & Kuhi 1979; Weaver & Jones 
1992; Kenyon & Hartmann 1995). Apart from the SED flux 
at 12 microns due to the 0.085Mj planet (dotted line) be- 
ing slightly higher than that observed, the SEDs computed 
from the density profiles fit the 2, 3, and 12 micron data 
points well. The flux between 3 and 12 microns varies quite 
considerably for the different masses, suggesting that more 
observations in this wavelength range could allow for the de- 
termination of the companion mass or, since the Reynolds 
number in the disc is unlikely to be accurately known, the 
relationship betwen the mass ratio and the Reynolds num- 
ber. The SEDs due to the 21 and 43Mj substellar objects 
also have a flux at 25 microns that is considerably higher 
than that observed. Although not definitive, this may sug- 
gest that if the GM Aur SED is due to an embedded com- 
panion (planet) orbiting at ~ 2.5 au, the mass must be such 
as to produce a gap that is not strongly truncated at the 
outer Lindblad resonance, i.e. q > TiT 1 but q 2 < 1ZT 1 . This 
also appears to be the case for TW Hydra as Calvet et al. 
(2002) require low density material within 4 au to fit the 
2-10 micron SED flux. 

At the distance of GM Aur (140 pc), a few Jupiter mass 
planet orbiting at 2.5 au will produce a 'wobble' in the cen- 
tral star with an angular displacement of ~ 0.1 milliarcsec, 
and a period of 4 years. Future high precision astrometric 
instruments should be sufficiently sensitive to make such an 
observation. Observations of this kind have been suggested 
for detecting both the presence of embedded planets (Boss 
1998) and for detecting instabilities in self-gravitating pro- 
toplanetary discs (Rice et al. 2002). 

Although it is difficult in SPH both to determine and 
set the viscosity, by performing a number of simulations with 
various companion masses we are able to determine how the 
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Figure 4. GM Aur's SED computed using the azimuthally aver- 
aged density profiles due to planets with masses of 0.085 (dotted 
line), 1.7 (solid line) and substellar companions with masses of 21 
(dash-dot line) and 43Mj (dashed line), as input to the Monte 
Carlo radiation transfer. 



radius (au) 

Figure 6. Azimuthally averaged midplane density profiles for a 
disc modelled, between 0.25 and 300 au using 300000 SPH par- 
ticles (solid line), and one modelled between 0.25 and 25 au us- 
ing 150000 SPH particles. The higher resolution simulation has 
a larger Reynolds number and hence, even though both simu- 
lations have the same planet mass (1.7 Mj), the gap is slightly 
harder than the lower resolution simulation. 




10 ' , Lj , I 

1 10 100 

Figure 5. Detailed view of the computed SEDs between 5 and 
200 microns. The line styles correspond to embedded planet 
masses of 0.085 (dotted line), 1.7 (solid line), 21 (dash-dot line), 
and 43Mj (dashed line). Error bars have been added to the data 
points with the error (±25 %) based on worst case error estimates. 

density profile, and resulting SED, varies for various val- 
ues of the mass ratio to Reynolds numbers ratio. To further 
test the relationship between planet mass and viscosity we 
have also performed a high resolution SPH simulation us- 
ing 150000 SPH particles to simulate the inner 25 au of 
the accretion disc. This increases the Reynolds number by 
a factor of 2 and allows the inner few au of the disc to be 
represented by more SPH particles once the gap has formed. 
Figure 6 shows the midplane disc density for the simulation 



using 300000 particles to represent the entire disc (solid line) 
and the simulation using 150000 particles to respresent the 
inner 25 au (dashed line). In both cases a planet mass of 
l.lMj was used. Although the density profiles differ, this is 
entirely consistent with the increased Reynolds number of 
the high resolution simulation. Extrapolating the high reso- 
lution simulation to 300 au and using the resulting density 
profile as input to the radiation transfer produces an SED 
that, like the previous 0.085 and 1.7Mj simulations, fits the 
data remarkably well. 



4 SUMMARY 

We have presented an SPH simulation of a planet opening up 
a gap in GM Aur's protoplanetary disk. For a 1.7Mj planet 
orbiting at 2.5 au, the disc density within 4 au is decreased 
by three orders of magnitude. The almost Keplerian velocity 
structure matches that derived from CO observations. Fur- 
ther, our radiation transfer simulations of a passively heated 
disk show that the disc structure from the SPH simulation 
reproduces the observed SED of GM Aur. 

Additional simulations of companions (planets) with 
masses of 0.085, 21, and 43Mj were also performed. The 
structure of the resulting gaps depend on the relationship 
between the mass ratio and the Reynolds number in the 
disc. For the Reynolds number in our simulations, the two 
planets (0.085 and 1.7Mj) produce gaps with relatively soft 
edges, while the two substellar objects (21 and 43Mj) trun- 
cate the disc strongly at the outer Lindblad resonance (~ 4 
au). As with the 1.7Mj planet, the SED computed using 
the density profile due to the 0.085Mj planet fits the data 
points extremely well. The SEDs computed using the den- 
sity profiles due to the two substellar objects (21 and 43Mj) 
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fit most data points well but have quite a significant excess 
at 25 microns. Although not definitive, this may suggest 
that the mass must be such as to produce a gap that is not 
strongly truncated at the outer Lindblad resonance. Further 
observations between 3 and 20 microns, where there is a con- 
siderable difference in the computed SEDs, but currently no 
data points, will be needed to determine if this is indeed the 
case. Future SIRTF observations will provide the required 
measurements to test our predictions. 

Since the density profile depends on the relationship 
between the mass ratio and the Reynolds number, addi- 
tional information would be required to determine the actual 
planet mass. A few Jupiter mass planet orbiting at 2.5 au 
will produce a 'wobble' in the central star, at the distance 
of GM Aur (140 pc), with an angular displacement of ~ 0.1 
milliarcsec, and a period of 4 years. Future high precision 
astrometric instruments should be sufficiently sensitive to 
make such observations. Not only would this establish if a 
planet is indeed orbiting at ~ 2.5 au but should also allow 
for the determination of the planet mass from which the 
Reynolds number, and hence the turbulent viscosity, could 
be estimated. 

The transfer of angular momentum between the em- 
bedded planet and the disc material is ultimately expected 
to cause the planet to migrate inwards (Lin & Paploizou 
1986). Simulations (R. Nelson et al 2000) suggest that the 
migration timescale is of order 10 5 yrs for a Jupiter mass 
planet starting at 5 au. Since as many as 15 percent of Solar 
type stars are expected to have planets within 5 au, that 
there are only currently two T Tauri stars with SED sig- 
natures of close-in massive planets is somewhat surprising. 
This could reflect the relatively sparse observations in the 
relevant wavelength range (3 — 12 microns). Alternatively, 
fully evacuated inner holes may form in only a fraction of 
classical T Tauri stars with planets. An annular gap with a 
sizeable inner disc would lead to less obvious signatures in 
the SED. 



ACKNOWLEDGEMENTS 

We acknowledge financial support from a PPARC standard 
grant (WKMR) and Advanced Fellowship (KW). BW would 
like to acknowledge support from NASA's Long Term Space 
Astrophysics (LTSA) Research Program, NAG 5-8412. 

REFERENCES 

Adams F.C., Lada C.J., Shu F.H., 1987, ApJ, 312, 788. 
Armitagc P.J., Livio M., Lubow S.H., Pringle J.E., 2002, MNRAS, 
334, 248 

Artymowicz P., Lubow S.H., 1994, ApJ, 421, 651. 
Artymowicz P., Lubow S.H., 1996, ApJ, 467, L77. 
Bate M.R., Bonncll I.A., Bromm V., 2002, MNRAS, 336, 705. 
Bate M.R., Bonnell LA., Price N.M., 1995, MNRAS, 277, 362. 
Beckwith S.V.W., 1999, in "The Origin of Stars and Planetary 

Systems," ed. C.J. Lada & N.D. Kylafis, Kluwer, 579 
Beckwith S.V.W., Sargent A. I., 1991, ApJ, 381, 205 
Beckwith S.V.W., Sargent A.I., Chini R.S., Gusten R., 1990, AJ, 

99, 924 

Benz W., 1990, in Buchler J.R., ed, The Numerical Modeling of 
Nonlinear Stellar Pulsations, Kluwer, Dordrecht, p. 269. 



Binncy J., Tremaine S., 1987, Galactic Dynamics, Princeton Uni- 
versity Press, Princeton. 
Bjorkman J.E., Wood K., 2001, ApJ, 554, 615 
Boss A. P., 1998, Nature, 393, 141. 
Boss A. P., 2000, ApJ, 553, 174. 

Calvet N., D'Alessio P., Hartmann L., Wilner D., Walsh A., Sitko, 

M., 2002, ApJ, 568, 1008 
Cardclli J.A., Clayton G.C., Mathis J.S., 1989, ApJ, 345, 245 
Chiang E.I., Goldreich P., 1999, ApJ, 519, 279 
Clarke C.J., Gcndrin A., Sotomayor M., 2001, MNRAS, 328, 485. 
Cohen M., Kuhi L.V., 1979, ApJS, 41, 743 

D'Alessio P., Canto J., Calvet N., Lizano S., 1998, ApJ, 500, 411. 
D'Alessio P., Calvet N., Hartmann L., Lizano S., Canto, J., 1999, 
ApJ, 527, 893 

D'Angelo G., Henning T., Kley W., 2002, A&A, 385, 647 

Duquennoy A., Mayor M., 1991, A&A, 248, 485. 

Dutrcy A., Guilloteau S., Prato L., Simon M., Duvert G., Schuster 

K., Menard F., 1998, A&A, 338, L63. 
Ghosh P., Lamb F.K., 1978, ApJ, 233, L83. 
Goldreich P., Tremaine S., 1980, ApJ, 241, 425 
Gullbring E., Hartmann L., Briccno C, Calvet N., 1998, ApJ, 

492, 323. 

Halbwachs J.-L., Arcnou F., Mayor M., Udry S., Queloz D., 2000, 

A&A, 355, 581. 
Haisch K.E., Lada E.A., Lada C.J., 2001, ApJ, 553, L153 
Hartmann L., 2002, ApJ, in press. 
Kcnyon S.J., Hartmann L., 1987, ApJ, 323, 714 
Kenyon S.J., Hartmann L., 1995, ApJS, 101, 117. 
Koerner D.W., Sargent A.I., Beckwith S.V.W., 1993 Icarus, 106, 

2 

Kuchncr M.J., Lccar M., 2002, ApJ, 574, L87. 

Kurucz R.L., 1994, CD-ROM 19, Solar Model Abundance Model 

Atmospheres (Cambridge: SAO) 
Lin D.N.C., Papaloizou J., 1979a, MNRAS, 186, 799. 
Lin D.N.C., Papaloizou J, 1979b, MNRAS, 188, 191. 
Lin D.N.C., Papalizou J., 1986, ApJ, 309, 846. 
Lin D.N.C., Bodenheimer P., Richardson D.C., 1996, Nature, 380, 

606 

Lissaucr J. J., 1993, ARA&A, 31, 129. 

Lubow S.H., Sicbcrt M., Artymowicz P., 1999, ApJ, 526, 1001. 
Lyndcn-Bell D., Pringle J.E., 1974, MNRAS, 168, 603. 
Marcy G.W., Butler R.P., 2000, PASP, 112, 137. 
Mathis J.S., Rumpl W., Nordsicck K.H., 1977, apj, 217, 425. 
Matsuyama I., Johnstone D., Hartmann L., 2002, ApJ, in press. 
Mayor M., Queloz D., 1995, Nature, 378, 355. 
Monaghan J.J., 1992, ARA&A, 30, 543. 
Murray J.R., 1996, MNRAS, 279, 402. 
Navarro J.F., White S.D.M., 1993, MNRAS, 265, 271. 
Nelson A.F., Benz W., Ruzmaikina T.V., 2000, ApJ, 529, 357. 
Nelson R.P., Papaloizou J.C.B., Massct F., Kley W., 2000, MN- 
RAS, 318, 18 
Pringle J.E., 1981, ARA&A, 19, 137. 

Rice W.K.M., Armitage P.J., Bate M.R., Bonncll LA., 2003, MN- 
RAS, 338, 227. 

Schneider G., Wood K., Silverstonc M., Hines D.C., Koerner 
D.W., Whitney B., Bjorkman J.E., Lowrance P.J., 2002, AJ, 
submitted. 

Shakura N.J., Sunyacv R.A., 1973, A&A, 24, 337. 

Sicss L., Forcstini M., Bcrtout C, 1999, A&A, 342, 480. 

Simon M., Prato L., 1995, ApJ, 450, 824. 

Simon M., Dutrey A., Guilloteau S., 2000, ApJ, 545, 1034. 

Sycr D., Clarke C.J., 1996, MNRAS, 278, L23. 

Trilling D.E., Benz W., Guillot T., Lunine J.I., Hubbard W.B., 

Burrows A., 1998, ApJ, 500, 428 
Trilling D.E., Lunine J. I., Benz W., 2002, A&A, 394, 241. 
Weaver W.B., Jones G., 1992, ApJS, 78, 239. 
Wolf S., Gueth F., Henning T., Kley W., 2002, ApJ, 566, L97. 



© 0000 RAS, MNRAS 000, 000-000 



8 Rice et al. 

Wood K., Lada C.J., Bjorkman J.E., Kenyon S.J., Whitney B.A., 

Wolff M.J., 2002a, ApJ, 567, 1138 
Wood K., Wolff M.J., Bjorkman J.E., Whitney B.A., 2002b, ApJ, 

564, 887 



© 0000 RAS, MNRAS 000, 000-000 



